The dimensions of the stacked image are 7,811 rows, 7,861 columns, 59,996,291 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.
The dimensions of the stacked image are 340 rows, 346 columns, 117,640 cells and 6 layers. The CRS is UTM on the WGS84 datum. The resolution is 30m x 30m.
A we can see, the NWDI, MNDWI, and WRI are all very similar in the way they pick up water cells with the same color, but they deviate differently in the dry ground cells. NDVI picks up water in a different color, and land in the opposite as well. The SWI seems to pick up only water cells and no land cells.
SWI has multiple NA values so we will have to replace those.
This shoes us they are extracted as a matrix of size 117640 x 6
##
## 1 2 3 4 5 6
## 0 488 17474 8361 53172 22289 655
## 1 13387 0 50 108 1646 10
| Flood Area m^2 | |
|---|---|
| NDVI | 5999400 |
| NDWI | 6490800 |
| MNDWI | 10745100 |
| WRI | 7622100 |
| SWI | 13680900 |
| K.Means | 12487500 |
Some of the values in the raster are not whole numbers because the color palette is continuous, so to make the overall raster, which is discrete, fit this color scheme, it is possible it combines values nearby and makes a raster value a non-whole number, which results in the color smoothing of the map.
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
library(tidyverse)
library(sf)
library(raster)
library(getlandsat)
library(mapview)
library(osmdata)
library(knitr)
bb = read_csv('../data/uscities.csv') %>%
filter(city == 'Palo') %>%
st_as_sf(coords = c('lng','lat'), crs = 4326) %>%
st_transform(5070) %>%
st_buffer(5000) %>%
st_bbox() %>%
st_as_sfc()
meta = read_csv('../data/palo-flood.csv')
files = lsat_scene_files(meta$download_url) %>%
filter(grepl(paste0('B',1:6,'.TIF$', collapse = '|'), file)) %>%
arrange(file) %>%
pull(file)
st = sapply(files, lsat_image)
b = stack(st) %>% setNames(paste0('band', 1:6))
cropper = bb %>% st_as_sf() %>% st_transform(crs(b))
cr_b = crop(b, cropper)
r <- cr_b %>% setNames(c('coastal', 'blue', 'green', 'red', 'nir', 'swir'))
plotRGB(r, r = 4, g= 3, b = 2)
plotRGB(r, r = 5, g= 4, b = 2)
plotRGB(r, r = 5, g= 6, b = 4)
plotRGB(r, r = 7, g= 5, b = 3)
plotRGB(r, r = 4, g= 3, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 4, b = 2, stretch = 'hist')
plotRGB(r, r = 5, g= 6, b = 4, stretch = 'lin')
plotRGB(r, r = 7, g= 5, b = 3, stretch = 'lin')
ndvi = (r$nir - r$red)/(r$nir + r$red)
ndwi = (r$green - r$nir)/(r$green + r$nir)
mndwi = (r$green - r$swir)/(r$green + r$swir)
wri = (r$green + r$red)/(r$nir + r$swir)
swi = 1/(sqrt(r$blue-r$swir))
rs <- stack(ndvi, ndwi, mndwi, wri, swi) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(rs, col = colorRampPalette(c("blue", "white", "red"))(256))
ndvi_t = function(x) {ifelse(x<=0,1,0)}
ndwi_t = function(x) {ifelse(x>=0,1,0)}
mndwi_t = function(x) {ifelse(x>=0,1,0)}
wri_t = function(x) {ifelse(x>=1,1,0)}
swi_t = function(x) {ifelse(x<=5,1,0)}
ndvi_f = calc(ndvi, ndvi_t)
ndwi_f = calc(ndwi, ndwi_t)
mndwi_f = calc(mndwi, mndwi_t)
wri_f = calc(wri, wri_t)
swi_f = calc(swi, swi_t)
f_stack = stack(c(ndvi_f, ndwi_f, mndwi_f, wri_f, swi_f)) %>% setNames(c('NDVI', 'NDWI', 'MNDWI', 'WRI', 'SWI'))
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
#sum(is.na(values(ndvi_f)))
#sum(is.na(values(ndwi_f)))
#sum(is.na(values(mndwi_f)))
#sum(is.na(values(wri_f)))
#sum(is.na(values(swi_f)))
values(swi_f) <- ifelse(is.na(values(swi_f)), 0, values(swi_f) )
#sum(is.na(values(swi_f)))
vals <- getValues(r)
#dim(vals)
vals <- na.omit(vals)
vs <- scale(vals)
idx <- 1:dim(vals)[1]
set.seed(20200905)
rast.clust <- kmeans(vs, 6, iter.max = 100)
kmeans_raster <- rs$MNDWI
values(kmeans_raster) <- rast.clust$cluster
plot(kmeans_raster)
table(getValues(swi_f), getValues(kmeans_raster))
k_thresh = function(x) {ifelse(x==1,1,0)}
k_flood = calc(kmeans_raster, k_thresh) %>% setNames('K-Means')
f_stack <- addLayer(f_stack, k_flood)
plot(f_stack, col = colorRampPalette(c("white","blue"))(256))
flood_area <- cellStats(f_stack, 'sum') * (30^2)
flood_area_km <- flood_area*0.000001
#names(flood_area)
kable(data.frame(flood_area), col.names = 'Flood Area m^2', caption = 'Flood Area for Different Methods')
rast_uncert <- calc(f_stack, function(x){sum(x)}) %>% setNames('Flood Uncertainty')
plot(rast_uncert, col = blues9, main = 'Flood Uncertainty')
values(rast_uncert) <- ifelse(values(rast_uncert)==0, NA, values(rast_uncert))
mapview(rast_uncert)